function [parms_fitted, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
    FitModel20220926_ver2(parms, kappa_guess, b_guess, eta_guess, z1_guess, chi_guess, ...
        U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, ave_yld_vol_target, NumTries, fsolve_opts)
% target moments:
% E[U], vol[u], vol(N*w)/vol(Y), vol(dlogY), vol(ylds), vol(mu3_1y)
% parameters:
% kappa, b, eta, z(1), chi, sigma_x_REC
% sigma_x(N,z) = exp(a + sigma_x_REC*1{z==1})

    if z1_guess>=0; error('z1_guess should be <0'); end
    if eta_guess<=0 || eta_guess>=1; error('check eta_guess'); end
    if kappa_guess<=0; error('check kappa0_guess'); end
    if b_guess<=0; error('check b_guess'); end
    if chi_guess>1 || chi_guess<0; error('check chi_guess'); end

    if parms.Nz~=2; error('Code is for 2 states.'); end
    if abs(parms.grid_z(2)-0)>1e-10; error('Code assumes z(2)=0.'); end
    
%     fixedpt_tol = 1e-8;
    % fixedpt_tol = 1e-7;
    fixedpt_tol = 1e-6; % speed things up
    
    N_Nxz = repmat(parms.grid_N(:),[1 parms.Nx parms.Nz]);
    N_Nz = repmat(parms.grid_N(:),[1 parms.Nz]);
    U_Nz = 1 - N_Nz;
    U_Nxz = 1 - N_Nxz;
    % x_Nxz = permute(repmat(parms.grid_x(:),[1 parms.NN parms.Nz]),[2 1 3]);
    
    % phi_of_x = @(x) 1./(1 + exp(-x));
    % phi_Nxz = phi_of_x(x_Nxz);
    % log_phi_Nxz = log(phi_Nxz);
    % input_mu3_Nxz = (1 - U_Nxz).*U_Nxz.*(1 - 2*U_Nxz).*(log_phi_Nxz.^3);
    
    FLAG_SUCCESSFUL_FSOLVE = false;
    
    for iter = 1:NumTries
        fprintf('iter %g with kappa_guess=%g, b_guess=%g, eta_guess=%g, z1_guess=%g, chi_guess=%g.\n', ...
            iter, kappa_guess, b_guess, eta_guess, z1_guess, chi_guess);
        [XFit, FitVal, EXITFLAG] = fsolve(@FitFun,[log(kappa_guess); ...
            log(b_guess); log(1/eta_guess-1); log(-z1_guess); ...
            log(1/chi_guess - 1)],fsolve_opts);
        if EXITFLAG>0
            fprintf('Found roots on iteration %g.\n', iter);
            FLAG_SUCCESSFUL_FSOLVE = true;
            break;
        else
            fprintf('iter %g, fitted vals (not solution):\n kappa=%g, b=%g, eta=%g, z1=%g, chi=%g.\n', ...
                iter, exp(XFit(1)), exp(XFit(2)), 1/(1 + exp(XFit(3))), ...
                -exp(XFit(4)), 1/(1 + exp(XFit(5))));
            kappa_guess = 0.1 + 5*rand;
            b_guess = 0.3 + 0.6*rand;
            eta_guess = 0.1*rand;
            z1_guess = fit_std_dz_2state(parms,vol_gy_target/sqrt(2.5+1*rand));
            chi_guess = 1/(1+4*rand);
        end
    end
    
    % do a last compute
    parms_fitted = parms;
    parms_fitted.kappa(:) = exp(XFit(1));
    parms_fitted.b(:) = exp(XFit(2));
    parms_fitted.eta(:) = 1/(1 + exp(XFit(3)));
    parms_fitted.grid_z(1) = -exp(XFit(4));
    parms_fitted.chi = 1/(1 + exp(XFit(5)));

    fprintf('kappa=%g, b=%g, eta=%g, z1=%g, chi=%g.\n', ...
        parms_fitted.kappa(1), parms_fitted.b(1), parms_fitted.eta(1), ...
        parms_fitted.grid_z(1), parms_fitted.chi);

    gesol = SolveFixedPoint_Model_20220926(parms_fitted,fixedpt_tol);
    parms_fitted.FitVal = FitVal;

    if ~FLAG_SUCCESSFUL_FSOLVE
        warning('Model fit: not all moments are satisfied.');
    end

    function Z = FitFun(X)
        
        Z = NaN*zeros(5,1);
        
        parms_temp = parms;
        parms_temp.kappa(:) = exp(X(1));
        parms_temp.b(:) = exp(X(2));
        parms_temp.eta(:) = 1/(1 + exp(X(3)));
        parms_temp.grid_z(1) = -exp(X(4));
        parms_temp.chi = 1/(1 + exp(X(5))); 
        
        z_Nxz = permute(repmat(parms_temp.grid_z(:),[1 parms.NN parms.Nx]),[2 3 1]); % note: this must be updated when z changes!!
        Y_Nxz = exp(z_Nxz).*N_Nxz;

        gesol_temp = SolveFixedPoint_Model_20220926(parms_temp,fixedpt_tol);
        if gesol_temp.VFIerr>fixedpt_tol
            warning('FitModel20220926_ver1: VFI error = %g is greater than tolerance of %g.', ...
                gesol_temp.VFIerr, fixedpt_tol)
            bComputeMoments = false;
        else
            bComputeMoments = true;
        end
        
        %% Newer version
        
        if bComputeMoments

        [Phi_Nxz,ss_dist_Nxz] = MakeSparseTransMatrix_Nxz(parms_temp,gesol_temp);

        % NEW CODE: quarterly average
        [mean_U, std_U] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,U_Nxz);
        [mean_Y, std_Y] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,Y_Nxz);
        [mean_Nw, std_Nw] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,N_Nxz.*gesol_temp.wages);
        ratio_stdNw_stdY = std_Nw/std_Y;
        
        % quarterly output vol (exact), from library\exact moments
        vol_gy = moment_quarterly_output_growth_vol_Nxz(parms_temp,gesol_temp,Phi_Nxz,ss_dist_Nxz);
        
        gesol_temp = PriceRealBonds_Model20220926(parms_temp,gesol_temp);
        
        % real yield curve
        mats = 12:12:60;
        real_yld_mean = zeros(size(mats));
        real_yld_std = zeros(size(mats));
        for n = 1:numel(mats)
            idx = find(gesol_temp.real_bonds.maturities==mats(n));
            mat_yr = mats(n)/12;
            [real_yld_mean(n),real_yld_std(n)] = Calc_monthly_vol(ss_dist_Nxz,-log(gesol_temp.real_bonds.prices(:,:,:,idx))/mat_yr);
        end
        real_yld_std = 100*real_yld_std;
        
        % std_mu3_cons_1y = CalcTimeSeriesDiffVol(input_mu3_Nxz,ss_dist_Nxz,Phi_Nxz,12);
        
        Z(1) = (mean_U - U_mean_target)/U_mean_target;
        Z(2) = (std_U - U_std_target)/U_std_target;
        Z(3) = (ratio_stdNw_stdY - wagebillratio_target)/wagebillratio_target;
        Z(4) = (vol_gy - vol_gy_target)/vol_gy_target;
        Z(5) = (mean(real_yld_std) - ave_yld_vol_target)/ave_yld_vol_target;
        
        end
        
    end

end